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Report  for  AASERT  grant  N00014-91-J-4037 


AASERT  grant  N00014-91-J-4037  supported  the  research  of  two  graduate  stu¬ 
dents,  James  M.  Keiser  and  Robert  Cramer.  James  Reiser  has  completed  his  research 
and  defended  his  thesis  April  17,  1995.  The  results  of  his  research  were  described  in 
the  report  submitted  5/31/95.  At  about  the  same  time  (May,  1995)  a  one  year  no  cost 
extension  for  this  grant  has  been  approved  in  order  to  allow  Robert  Cramer  to  finish 
his  research  towards  Ph.D.  At  the  time  of  this  writing,  Robert  Cramer  has  achieved  the 
goals  of  his  thesis  work  and  is  in  the  process  of  finishing  his  thesis  (which  he  will  defend 
this  summer). 

The  work  of  Robert  Cramer  is  concerned  with  computing  particle  interactions 
in  a  manner  consistent  with  projection  methods.  Using  a  wavelet-based  approach,  it 
appears  possible  to  construct  a  high  order  generalization  of  local  correction  methods 
(local  correction  methods  typically  work  well  if  low  accuracy  is  required).  A  brief  outline 
is  presented  below. 

Summary  of  Robert  Cramer’s  thesis 

The  work  of  Robert  Cramer  concerns  a  fast  wavelet  based  method  for  computing 
“particle  interactions”.  By  “particle  interactions”  we  understand  discrete  sums  repre¬ 
senting  interactions  between  points  in  space.  We  assume  that  the  interaction  between 
any  two  particles  depends  only  on  the  distance  between  them. 

We  wish  to  compute  the  sums 

N 

9{^m)  =  K{Xm  -  yn)f{yn)  ,  I  <  TU  <  M 

n=l 

Vn^^m 

in  a  number  of  operations  proportional  to  (M  -t-  N)  (rather  than  M  ■  N)  with  a  fixed  but 
arbitrary  accuracy.  In  the  context  of  a  particle  simulation,  the  function  f{y)  represents 
the  charge  density  sampled  at  the  discrete  set  of  points  {yi, . .  • ,  and  g{x)  represents 
a  potential  field  sampled  at  the  discrete  set  of  points  {xi, . . .  ,0:^/}.  The  requirement 
7^  yn  excludes  the  self-interaction  (which  is  generally  infinite). 
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It  is  often  possible  to  view  the  summation  problem  as  a  discretization  of  a  bounded 
integral  operator, 

g(x)=Kf(x)  =  jK(x-y)f(y)dy. 

However,  our  algorithms  do  not  depend  on  this  and,  moreover,  we  are  especially  interested 
in  kernels  for  which  this  interpretation  is  not  valid. 

The  applications  of  such  algorithms  are  numerous,  and  in  the  physical  sciences 
include  the  study  of  stellar  and  galactic  clusters,  gas  dynamics,  fluid  dynamics  (vortex 
methods),  flow  of  electrons  in  semi-conductor  devices,  and  the  study  of  ionized  liquids 
and  phase  changes  in  chemistry.  In  addition,  such  algorithms  permit  the  design  of  fast 
Fourier  transforms  for  unequally  spaced  data. 

A  fast  algorithm  to  compute  these  sums  is  provided  by  the  Fast  Multipole  Method 
(FMM).  We  are  looking  for  an  alternative  to  FMM  which  will  be  more  generic  and  fit 
naturally  into  projection-type  methods.  It  appears  that  the  wavelet  approach  leads  to  an 
algorithm  which  may  be  thought  of  as  a  high  order  generalization  of  the  method  of  local 
corrections.  These  local  corrections  methods  have  largely  been  of  relatively  low  order 
in  accuracy  (roughly  10“^  to  and  have  been  restricted  to  kernels  arising  from 

diflferential  equations,  mainly  Laplace’s  or  Poisson’s.  Our  approach  yields  methods  of 
much  broader  scope,  including  arbitrary  order  accuracy,  and  the  ability  to  handle  a  wider 
range  of  kernels.  Using  a  multiresolution  approach  it  appears  possible  to  analytically 
construct  an  efficient  high  order  approximation  of  the  correction  operator.  The  approach 
is  efficient  and  should  be  competetive  with  FMM. 

The  fact  which  provides  a  common  foundation  for  fast  summation  algorithms  is 
that  kernels  of  the  type  mentioned  above  can  be  split  into  the  sum  of  a  smooth  kernel 
which  describes  the  long-range  interactions,  and  a  singular  kernel  which  influences  only 
the  short-range  interactions.  Thus,  the  strategy  is  to  construct  a  globally  smooth  kernel, 
which  approximates  K{x  —  y)  well  if  the  distance  |a;  —  y\  is  large,  and  then  define 
the  short-range  part  through  the  equation 

^(SR)  —  _  j^(LR) 

The  long-range  part,  being  globally  smooth  nnd  a  convolution,  can  be  efficiently  applied 


using  the  Fourier  transform. 

The  singular  or  “high-frequency”  part  of  the  kernel  is  applied  directly,  but  due  to 
the  rapid  decay  of  in  the  physical  space,  these  short-range  interactions  involve  only 
a  few  particles.  For  this  reason,  this  stage  of  the  computation  can  also  be  accomplished 
at  relatively  low  computational  cost. 

Let  us  first  describe  the  one-dimensional  version  of  our  algorithm  and  then  indicate 
the  modifications  necessary  to  extend  this  approach  to  two  or  three  dimensions.  Let  us 
consider  a  multiresolution  analysis,  i.e.,  a  chain  of  closed  subspaces  of  L^(R), 

•  •  •  C  V2  C  Vi  C  Vo  C  K.1  C  VL2  C  •  •  • . 

The  approximation  properties  of  a  MRA  follow  from  two  essential  ingredients.  First, 
the  basis  functions  ^  are  compactly  supported,  and  this  support  becomes  arbitrarily 
small  as  we  go  to  finer  and  finer  subspaces.  Second,  polynomials  up  to  a  given  degree 
can  be  represented  exactly  in  the  MRA.  Namely,  with  each  MRA.  there  is  an  associated 
parameter  M,  and  for  each  m  =  0, . . . ,M,  there  exist  coefficients  G  R,  A;  G  Z  such 
that 

=  E44W- 

k 

Using  the  properties  of  compact  support  and  the  ability  to  accurately  represent 
functions  which  behave  like  polynomials  on  compact  intervals,  it  is  possible  to  construct 
a  kernel  Tj  which  operates  on  the  subspace  Vj,  and  which  approximates  K  as  follows. 
Given  any  e  >  0,  there  exists  a  constant  B  =  B{e),  and  an  integer  j  <  0,  such  that 

\K{x  —  y)—  Tj{x,  y)\  <  e  whenever  \x  —  y\>  2^B  . 

Our  multiresolution  kernel  has  the  explicit  form 


k  I 


where  the  coefficients  are  given  by 

^k-i=  J  J  K{x-  y)(tPk{x)4^i  (y)  dydx . 


The  short-range  part  is  then  defined  by 

K{x-y)-Tj{x,y)  ,\x-y\<2iB 

0  ,  k  -  y|  >  . 

Computing  the  long-range  contribution  to  the  potential  is  reduced  to  the  problem  of 
applying  a  Toeplitz  matrix  to  a  vector  which  is  accomplished  using  the  FFT.  To  evaluate 
the  short-range  contribution,  it  is  necessary  to  evaluate  the  kernel  Tj{x,y)  at  a  given 
point  {x,y).  For  this  we  provide  a  Fourier  series  expansion  of  the  kernel, 

y)  =  £  cos  |^2n7r  (^^)]  -  y) . 

This  series  converges  very  quickly,  and  we  typically  retain  only  the  first  four  terms.  The 
functions  are  expressed  in  terms  of  the  basis  function  and  are  tabulated 

for  rapid  evaluation.  Application  of  the  short-range  part  of  the  kernel  is  thus  an  0{M) 
procedure,  as  for  a  given  point  Xm,  I  <  m  <  M ,  only  a  few  of  the  N  particles  will  be 
involved. 

In  two  (and  higher)  spatial  dimensions,  we  approximate  the  kernel  K{x-x',  y-xf) 
by  constructing  an  MRA  kernel 


Tj{x,y,x',y')  =  J2i-k',i-i'4>{x')4{y') 

ic,i  y,v 

to  satisfy 

\K{x-x\y-y')-Tj{x,y,x\^)\<e 

whenever 

mQx{\x-x'\,\y-y'\}>2^B 
for  appropriately  chosen  B  =  B{e)  and  j  <  0. 

The  long-range  contribution  is  computed  using  the  FFT,  entirely  analogous  to 
the  one-dimensional  case.  An  additional  step  required  in  higher  dimensions  in  order 
to  rapidly  evaluate  Tj{x,y,x\y')  at  a  point  {x,y,x',y'),  is  the  step  of  computing  the 
singular-value  decomposition  of  the  coefficient  matrix  to  obtain  the  representation 


where  {a^}  are  the  singular  values  and  R  is  the  numerical  rank.  Due  to  the  fact  that  the 
coefficient  matrix  is  of  low  rank,  the  parameter  R  is  relatively  small  (for  example, 
i?  =  8  for  single  precision).  We  obtain  the  two-dimensional  equation  for  evaluating  the 
short-range  interactions. 


The  functions  and  are  again  tabulated  for  rapid  evaluation.  The  number  of 
terms  necessary  to  retain  in  the  sums  above  is  very  small  (typically  four)  since  the  series 
converges  rapidly. 

Further  details  may  be  obtained  upon  request  from  Robert  Cramer, 
bcramer@boulder.colorado.edu. 
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